#this program demonstrates some of the network randomization bits and
#example models.
#clear everything to start...
rm(list = ls())
gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 375352 20.1 750400 40.1 592000 31.7
## Vcells 583261 4.5 1308461 10.0 786378 6.0
#load basic data manipulation bits
library(dplyr);
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(readr);
library(magrittr) #Supports pipe (%>%) commands that allow you to perform multiple operations with one statement
library(tidyr) #Additional functions for manipulating data
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
library(ggplot2)
#load the data - this is the same code as before
setwd("/Users/mariacristinaramos/Documents/MARIA CRISTINA/ACTIVE/SN-H2018/")
AHS_Base <- read_csv ("ahs_wpvar (1).csv",
col_names = TRUE);
## Parsed with column specification:
## cols(
## .default = col_integer(),
## POP_BEH = col_double(),
## PSICK_SS = col_double(),
## PSICK_H = col_double(),
## PSICK_L = col_double(),
## POP_BEH_SS = col_double()
## )
## See spec(...) for full column specifications.
AHS_adjlist <- AHS_Base %>%
select(ego_nid, mfnid_1:mfnid_5, ffnid_1:ffnid_5, grade, sex, PSMOKES,commcnt) %>%
filter(commcnt==1);
AHS_Edges <- AHS_adjlist %>%
rename( id = `ego_nid`,
gender = `sex`) %>%
gather(Alter_Label, Target, mfnid_1:mfnid_5, ffnid_1:ffnid_5, na.rm = TRUE)
AHS_Edges=AHS_Edges %>% filter (Target != 99999);
AHS_Edges=AHS_Edges %>%select(id, Target);
#install.packages("statnet")
library(statnet)
## Loading required package: tergm
## Loading required package: statnet.common
##
## Attaching package: 'statnet.common'
## The following object is masked from 'package:base':
##
## order
## Loading required package: ergm
## Loading required package: network
## network: Classes for Relational Data
## Version 1.13.0 created on 2015-08-31.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Martina Morris, University of Washington
## Skye Bender-deMoll, University of Washington
## For citation information, type citation("network").
## Type help("network-package") to get started.
##
## ergm: version 3.8.0, created on 2017-08-18
## Copyright (c) 2017, Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Carter T. Butts, University of California -- Irvine
## Steven M. Goodreau, University of Washington
## Pavel N. Krivitsky, University of Wollongong
## Martina Morris, University of Washington
## with contributions from
## Li Wang
## Kirk Li, University of Washington
## Skye Bender-deMoll, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").
## NOTE: Versions before 3.6.1 had a bug in the implementation of the
## bd() constriant which distorted the sampled distribution somewhat.
## In addition, Sampson's Monks datasets had mislabeled vertices. See
## the NEWS and the documentation for more details.
## Loading required package: networkDynamic
##
## networkDynamic: version 0.9.0, created on 2016-01-12
## Copyright (c) 2016, Carter T. Butts, University of California -- Irvine
## Ayn Leslie-Cook, University of Washington
## Pavel N. Krivitsky, University of Wollongong
## Skye Bender-deMoll, University of Washington
## with contributions from
## Zack Almquist, University of California -- Irvine
## David R. Hunter, Penn State University
## Li Wang
## Kirk Li, University of Washington
## Steven M. Goodreau, University of Washington
## Jeffrey Horner
## Martina Morris, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("networkDynamic").
##
## tergm: version 3.4.1, created on 2017-09-12
## Copyright (c) 2017, Pavel N. Krivitsky, University of Wollongong
## Mark S. Handcock, University of California -- Los Angeles
## with contributions from
## David R. Hunter, Penn State University
## Steven M. Goodreau, University of Washington
## Martina Morris, University of Washington
## Nicole Bohme Carnegie, New York University
## Carter T. Butts, University of California -- Irvine
## Ayn Leslie-Cook, University of Washington
## Skye Bender-deMoll
## Li Wang
## Kirk Li, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("tergm").
## Loading required package: ergm.count
##
## ergm.count: version 3.2.2, created on 2016-03-29
## Copyright (c) 2016, Pavel N. Krivitsky, University of Wollongong
## with contributions from
## Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm.count").
## NOTE: The form of the term 'CMP' has been changed in version 3.2
## of 'ergm.count'. See the news or help('CMP') for more information.
## Loading required package: sna
## sna: Tools for Social Network Analysis
## Version 2.4 created on 2016-07-23.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## For citation information, type citation("sna").
## Type help(package="sna") to get started.
##
## statnet: version 2016.9, created on 2016-08-29
## Copyright (c) 2016, Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Carter T. Butts, University of California -- Irvine
## Steven M. Goodreau, University of Washington
## Pavel N. Krivitsky, University of Wollongong
## Skye Bender-deMoll
## Martina Morris, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("statnet").
## unable to reach CRAN
g=as.network(AHS_Edges)
g %v% "grade" <- AHS_adjlist$grade
g %v% "sex" <- AHS_adjlist$sex
g %v% "smokes" <- AHS_adjlist$PSMOKES
# get degree scores
outdegree <- degree(g, cmode = "outdegree")
indegree<- degree(g, cmode = "indegree")
g %v% "indegree" <- indegree
g %v% "outdegree" <- outdegree
g %v% "degree" <- degree(g)
plot(g, vertex.col="grade")
plot(g, vertex.col="sex")
#how transitive is our network?
obs_tran=gtrans(g)
obs_tran
## [1] 0.375
#distribution of dyads
dyads=dyad.census(g)
dyads
## Mut Asym Null
## [1,] 85 135 2265
#to simulate...first initialize an empty vector
simsize=500;
ran_tran=vector(length=simsize)
#loop for the simulation...
j=1
while(j<(simsize+1)){
r=rguman(1,nv=71,mut=dyads[1],asym=dyads[2],null=dyads[3])
ran_tran[j]=gtrans(r)
if (is.na(ran_tran[j])) next else j=j+1 #skip bad sim cases
}
dt=as.data.frame(ran_tran)
ggplot(dt,aes(x=ran_tran))+geom_histogram()+
geom_vline(xintercept=obs_tran,
linetype=1,
color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#perhaps the clustering observed is due to general categorical mixing,
#by grade or sex?
grades=as.character(g %v% "grade")
grdmix=mixingmatrix(g, "grade")
grdmix
## To
## From 7 8 9 10 11 12 Total
## 7 52 5 1 1 0 0 59
## 8 8 33 9 0 1 1 52
## 9 0 10 70 1 4 1 86
## 10 0 0 3 30 10 0 43
## 11 1 0 2 7 43 4 57
## 12 0 0 1 0 2 5 8
## Total 61 48 86 39 60 11 305
#should be able to do this with the RGNMIX function, but that's
#not working, so let's do the same thing by hand.
#here's my failed code for the build in function,
################
#to simulate...first initialize an empty vector
#simsize=500;
#ran_tran=vector(length=simsize)
#loop for the simulation...
#j=1
#while(j<(simsize+1)){
# ran_g<-rgnmix(1,grades,grdmix$matrix,method="exact",return.as.edgelist=TRUE)
# ran_tran[j]=gtrans(as.network(ran_g))
# j=j+1
# if (is.na(ran_tran[j])) next else j=j+1 #skip bad sim cases
#}
#dt=as.data.frame(ran_tran)
#ggplot(dt,aes(x=ran_tran))+geom_histogram()+
# geom_vline(xintercept=obs_tran,
# linetype=1,
# color="red")
#now do it by hand. This code steals directly from some stuff
#jake fisher put together in 2017 session.
#1. get a probability matrix from the mixing matrix...need a denominator
c=as.matrix(table(grades))
cells=c%*%t(c)
mixprob=grdmix$matrix/cells
mixprob
## To
## From 7 8 9 10 11
## 7 0.520000000 0.038461538 0.025000000 0.006666667 0.000000000
## 8 0.061538462 0.195266272 0.173076923 0.000000000 0.005917160
## 9 0.000000000 0.192307692 4.375000000 0.016666667 0.076923077
## 10 0.000000000 0.000000000 0.050000000 0.133333333 0.051282051
## 11 0.007692308 0.000000000 0.038461538 0.035897436 0.254437870
## 12 0.000000000 0.000000000 0.015625000 0.000000000 0.009615385
## To
## From 12
## 7 0.000000000
## 8 0.004807692
## 9 0.015625000
## 10 0.000000000
## 11 0.019230769
## 12 0.019531250
#note this is technically not quite right...as we should subtractk
#1 from the diagonal...but I'm being lazy...
# 2. Now create an empty matrix:
prob.matrix <- matrix(NA, ncol = network.size(g),
nrow = network.size(g))
# 3. ...fill in the values based on the mixing matrix. Going to
#so get the categories...
grade.values <- sort(unique(grades))
# The general idea is to do this...
#prob.matrix[which(grades == 7), which(grades == 7)] <- tie.probs[1, 1]
# ...but that's going to require 36 lines of code. Let's write a loop instead.
for (i in 1:nrow(mixprob)) {
for (j in 1:ncol(mixprob)) {
prob.matrix[which(grades == grade.values[i]),
which(grades == grade.values[j])] <- mixprob[i, j]
}
}
#now we have a matrix where each ij cell is the probabilty of a
#tie based on mixing...so just do random bernulli draws from that,
#using the same code strucure as above...
simsize=500;
ran_tran=vector(length=simsize)
#loop for the simulation...
j=1
while(j<(simsize+1)){
r=rgraph(n = network.size(g), tprob = prob.matrix)
ran_tran[j]=gtrans(r)
if (is.na(ran_tran[j])) next else j=j+1 #skip bad sim cases
}
dt=as.data.frame(ran_tran)
ggplot(dt,aes(x=ran_tran))+geom_histogram()+
geom_vline(xintercept=obs_tran,
linetype=1,
color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#could do the same thing by Sex, but the plot we had suggested it's
#not strongly clustered by sex...but you get the idea...
#also, would be more efficient to create a sample of graphs, store them
#then you can test a bunch of stats agains the same set...but transitivyt
#is cheap to calculate...so using that...
plot(hist(ran_tran),xlim = c(0, .5))
abline(v = obs_tran, col = "red")
# Now we pass both of those to the function in igraph
#one quick test...
r=igraph::sample_degseq(out.deg = outdegree, in.deg = indegree,method = "simple.no.multiple")
igraph::transitivity(r)
## [1] 0.1429688
igraph::plot.igraph(r)
#same loop bit...could just write a function for this...
simsize=50;
ran_tran=vector(length=simsize)
#loop for the simulation...note I use the igraph transitivity function
j=1
while(j<(simsize+1)){
r=igraph::sample_degseq(out.deg = outdegree, in.deg = indegree,method = "simple.no.multiple")
ran_tran[j]=igraph::transitivity(r)
if (is.na(ran_tran[j])) next else j=j+1 #skip bad sim cases
}
dt=as.data.frame(ran_tran)
ggplot(dt,aes(x=ran_tran))+geom_histogram()+
geom_vline(xintercept=obs_tran,
linetype=1,
color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#again...not even close. But, all this is binary; let's build a model
#to look at this with multiple conditionals.
####################
# ERGM
####################
#ergm code adapted from Jeff Smith's tutorial.
#clean house...
rm(c,cells, dt, dyads, grdmix,prob.matrix,r)
rm(grade.values,grades,i,j,mixprob,obs_tran,rant_tran,simsize)
## Warning in rm(grade.values, grades, i, j, mixprob, obs_tran, rant_tran, :
## object 'rant_tran' not found
#super simple model..
mod.rand<-ergm(g~edges)
## Evaluating log-likelihood at the estimate.
##
summary(mod.rand)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges
##
## Iterations: 6 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -2.7275 0.0591 0 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 2293 on 4969 degrees of freedom
##
## AIC: 2295 BIC: 2302 (Smaller is better.)
#interpretation of model coefficients:
#the log-odds of any tie existing is -2.7275
#probability: exp(-2.7275)/(1+exp(-2.7275))=.0613
exp(-2.7275)/(1+exp(-2.7275))
## [1] 0.06137001
#which compares to density as..
gden(g)
## [1] 0.06136821
#How good is the model?
try(plot(gof(mod.rand), main=NULL))
#add some covariates...
mod.homoph<-ergm(g~edges+nodematch("sex")+nodematch("grade")
+nodematch('smokes'))
## Evaluating log-likelihood at the estimate.
summary(mod.homoph)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes")
##
## Iterations: 7 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -4.5787 0.1850 0 < 1e-04 ***
## nodematch.sex 0.4621 0.1311 0 0.000428 ***
## nodematch.grade 3.0493 0.1422 0 < 1e-04 ***
## nodematch.smokes 0.4062 0.1483 0 0.006196 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 1712 on 4966 degrees of freedom
##
## AIC: 1720 BIC: 1746 (Smaller is better.)
mod.gof=gof(mod.homoph)
plot(mod.gof, main=NULL)
#add some simple volume bitscovariates...
mod.p1<-ergm(g~edges+nodematch("sex")+nodematch("grade")
+nodematch('smokes')+sender+receiver)
## Observed statistic(s) sender4, sender5, sender15, sender17, sender42, sender57, sender63, receiver4, receiver17, and receiver44 are at their smallest attainable values. Their coefficients will be fixed at -Inf.
## Evaluating log-likelihood at the estimate.
summary(mod.p1)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") +
## sender + receiver
##
## Iterations: 7 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -6.54235 1.39273 0 <1e-04 ***
## nodematch.sex 0.72340 0.15509 0 <1e-04 ***
## nodematch.grade 3.60634 0.17189 0 <1e-04 ***
## nodematch.smokes 1.27256 0.20774 0 <1e-04 ***
## sender2 2.22577 1.25583 0 0.0764 .
## sender3 0.85358 1.31365 0 0.5159
## sender4 -Inf 0.00000 0 <1e-04 ***
## sender5 -Inf 0.00000 0 <1e-04 ***
## sender6 1.17264 1.29221 0 0.3642
## sender7 1.82007 1.28618 0 0.1571
## sender8 2.33419 1.27155 0 0.0665 .
## sender9 2.03254 1.27211 0 0.1102
## sender10 -0.12412 1.41704 0 0.9302
## sender11 -0.74302 1.60190 0 0.6428
## sender12 -1.01392 1.57899 0 0.5208
## sender13 1.44915 1.30637 0 0.2674
## sender14 1.33064 1.32471 0 0.3152
## sender15 -Inf 0.00000 0 <1e-04 ***
## sender16 2.34335 1.25453 0 0.0618 .
## sender17 -Inf 0.00000 0 <1e-04 ***
## sender18 1.99117 1.29061 0 0.1229
## sender19 -0.07322 1.40150 0 0.9583
## sender20 1.52207 1.28439 0 0.2361
## sender21 2.09979 1.26439 0 0.0968 .
## sender22 0.88576 1.34555 0 0.5104
## sender23 1.46870 1.28751 0 0.2540
## sender24 0.88576 1.34555 0 0.5104
## sender25 1.25861 1.34235 0 0.3485
## sender26 1.34868 1.27622 0 0.2907
## sender27 0.51817 1.35188 0 0.7015
## sender28 0.30073 1.35208 0 0.8240
## sender29 2.43662 1.28473 0 0.0579 .
## sender30 1.36028 1.32210 0 0.3036
## sender31 1.88595 1.28043 0 0.1408
## sender32 0.72578 1.32209 0 0.5831
## sender33 2.09444 1.25897 0 0.0963 .
## sender34 2.92916 1.26585 0 0.0207 *
## sender35 1.48498 1.28646 0 0.2484
## sender36 2.04011 1.27129 0 0.1086
## sender37 1.05055 1.32819 0 0.4290
## sender38 2.06804 1.28481 0 0.1075
## sender39 1.03670 1.40765 0 0.4615
## sender40 1.71434 1.26586 0 0.1757
## sender41 -0.86750 1.56951 0 0.5805
## sender42 -Inf 0.00000 0 <1e-04 ***
## sender43 1.97972 1.28468 0 0.1234
## sender44 1.84499 1.28992 0 0.1527
## sender45 0.37849 1.35551 0 0.7801
## sender46 1.20985 1.30077 0 0.3524
## sender47 -0.07322 1.40150 0 0.9583
## sender48 0.83024 1.30801 0 0.5256
## sender49 1.82994 1.29996 0 0.1593
## sender50 0.39527 1.35435 0 0.7704
## sender51 2.59763 1.28047 0 0.0425 *
## sender52 2.39671 1.28765 0 0.0628 .
## sender53 0.74770 1.30858 0 0.5678
## sender54 0.63491 1.34110 0 0.6359
## sender55 0.87978 1.41662 0 0.5346
## sender56 1.10414 1.30234 0 0.3966
## sender57 -Inf 0.00000 0 <1e-04 ***
## sender58 2.32900 1.29162 0 0.0714 .
## sender59 0.96191 1.31592 0 0.4648
## sender60 1.20635 1.29933 0 0.3532
## sender61 2.62317 1.27827 0 0.0402 *
## sender62 -0.79590 1.57732 0 0.6139
## sender63 -Inf 0.00000 0 <1e-04 ***
## sender64 1.37381 1.28712 0 0.2859
## sender65 2.34097 1.29990 0 0.0718 .
## sender66 2.35099 1.28906 0 0.0682 .
## sender67 0.63782 1.41165 0 0.6514
## sender68 0.77762 1.32296 0 0.5567
## sender69 1.78102 1.27294 0 0.1618
## sender70 0.09125 1.43162 0 0.9492
## sender71 0.86733 1.45067 0 0.5499
## receiver2 -1.43489 1.07362 0 0.1814
## receiver3 -1.69160 1.06305 0 0.1116
## receiver4 -Inf 0.00000 0 <1e-04 ***
## receiver5 -0.79045 0.97874 0 0.4193
## receiver6 -1.66543 1.06255 0 0.1171
## receiver7 -0.66723 0.95563 0 0.4851
## receiver8 0.22141 0.91166 0 0.8081
## receiver9 0.36172 0.86484 0 0.6758
## receiver10 -0.87903 0.95032 0 0.3550
## receiver11 -1.25341 1.00146 0 0.2108
## receiver12 -2.50486 1.27755 0 0.0500 *
## receiver13 -1.65921 1.08480 0 0.1262
## receiver14 -0.19915 0.93421 0 0.8312
## receiver15 -0.21006 0.90663 0 0.8168
## receiver16 -1.03081 0.99540 0 0.3005
## receiver17 -Inf 0.00000 0 <1e-04 ***
## receiver18 -1.00664 1.02771 0 0.3274
## receiver19 -0.94528 0.92624 0 0.3075
## receiver20 0.13942 0.88488 0 0.8748
## receiver21 -1.59147 1.08360 0 0.1420
## receiver22 -1.73458 1.09012 0 0.1116
## receiver23 -0.12357 0.88993 0 0.8896
## receiver24 -1.73458 1.09012 0 0.1116
## receiver25 -0.42157 0.98180 0 0.6677
## receiver26 -1.22260 0.96921 0 0.2072
## receiver27 0.32993 0.87520 0 0.7062
## receiver28 -2.47186 1.27910 0 0.0534 .
## receiver29 0.47522 0.91342 0 0.6029
## receiver30 0.12959 0.91004 0 0.8868
## receiver31 0.51347 0.87020 0 0.5552
## receiver32 -1.66313 1.06639 0 0.1189
## receiver33 0.31993 0.87070 0 0.7133
## receiver34 1.66643 0.85103 0 0.0503 .
## receiver35 0.12686 0.87443 0 0.8847
## receiver36 0.57683 0.85120 0 0.4980
## receiver37 -0.54970 0.94518 0 0.5609
## receiver38 0.20867 0.91396 0 0.8194
## receiver39 -0.30300 0.98849 0 0.7592
## receiver40 -1.12114 0.97563 0 0.2506
## receiver41 -1.32436 0.96960 0 0.1720
## receiver42 0.75364 0.87979 0 0.3917
## receiver43 -0.14799 0.92850 0 0.8734
## receiver44 -Inf 0.00000 0 <1e-04 ***
## receiver45 -1.25725 0.99393 0 0.2060
## receiver46 -0.15188 0.90265 0 0.8664
## receiver47 -0.94528 0.92624 0 0.3075
## receiver48 -0.34516 0.87911 0 0.6946
## receiver49 -0.28730 0.97578 0 0.7684
## receiver50 -0.47568 0.91303 0 0.6024
## receiver51 1.82942 0.84156 0 0.0298 *
## receiver52 0.59043 0.90601 0 0.5146
## receiver53 -1.74439 1.05338 0 0.0978 .
## receiver54 0.41471 0.84789 0 0.6248
## receiver55 -1.61996 1.28711 0 0.2082
## receiver56 -1.72440 1.07485 0 0.1087
## receiver57 -1.70193 1.08087 0 0.1154
## receiver58 0.05836 0.94460 0 0.9507
## receiver59 -0.63772 0.95194 0 0.5029
## receiver60 0.34174 0.86306 0 0.6921
## receiver61 1.19676 0.87435 0 0.1711
## receiver62 -1.68211 1.08277 0 0.1204
## receiver63 -0.43936 1.06709 0 0.6806
## receiver64 -2.49890 1.28600 0 0.0521 .
## receiver65 -0.37488 1.01075 0 0.7107
## receiver66 0.68781 0.90364 0 0.4466
## receiver67 -0.02045 0.94255 0 0.9827
## receiver68 -0.77370 0.94078 0 0.4109
## receiver69 0.15329 0.88341 0 0.8622
## receiver70 -1.75929 1.08335 0 0.1045
## receiver71 -0.47914 1.11343 0 0.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 2292 on 4826 degrees of freedom
##
## AIC: 2580 BIC: 3517 (Smaller is better.)
##
## Warning: The following terms have infinite coefficient estimates:
## sender4 sender5 sender15 sender17 sender42 sender57 sender63 receiver4 receiver17 receiver44
mod.gof=gof(mod.p1)
plot(mod.gof, main=NULL)
#model predicts better, but huge BIC cost! Simplify?
g %v% "indegree" <- indegree
g %v% "outdegree" <- outdegree
mod.p2<-ergm(g~edges+nodematch("sex")+nodematch("grade")
+nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
+mutual)
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20:
## Optimizing with step length 0.789958896639935.
## The log-likelihood improved by 3.242.
## Iteration 2 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.1611.
## Step length converged once. Increasing MCMC sample size.
## Iteration 3 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.02584.
## Step length converged twice. Stopping.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
summary(mod.p2)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") +
## nodeicov("indegree") + nodeocov("outdegree") + mutual
##
## Iterations: 3 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -8.33888 0.37137 0 <1e-04 ***
## nodematch.sex 0.53239 0.12932 0 <1e-04 ***
## nodematch.grade 2.60143 0.15791 0 <1e-04 ***
## nodematch.smokes 0.86950 0.15439 0 <1e-04 ***
## nodeicov.indegree 0.29062 0.02828 0 <1e-04 ***
## nodeocov.outdegree 0.32907 0.03489 0 <1e-04 ***
## mutual 2.32249 0.23160 0 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 1321 on 4963 degrees of freedom
##
## AIC: 1335 BIC: 1381 (Smaller is better.)
mod.gof=gof(mod.p2)
plot(mod.gof, main=NULL)
mod.gendmut<-ergm(g~edges+nodematch("sex")+nodematch("grade")
+nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
+mutual(same="sex",diff=TRUE))
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20:
## Optimizing with step length 0.953219017725786.
## The log-likelihood improved by 3.079.
## Iteration 2 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.3219.
## Step length converged once. Increasing MCMC sample size.
## Iteration 3 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.03451.
## Step length converged twice. Stopping.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
summary(mod.gendmut)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") +
## nodeicov("indegree") + nodeocov("outdegree") + mutual(same = "sex",
## diff = TRUE)
##
## Iterations: 3 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -8.18801 0.38824 0 <1e-04 ***
## nodematch.sex -0.14310 0.17442 0 0.412
## nodematch.grade 2.83441 0.15881 0 <1e-04 ***
## nodematch.smokes 0.95223 0.16382 0 <1e-04 ***
## nodeicov.indegree 0.30719 0.02848 0 <1e-04 ***
## nodeocov.outdegree 0.35012 0.03357 0 <1e-04 ***
## mutual.same.sex.1 2.77943 0.37059 0 <1e-04 ***
## mutual.same.sex.2 2.53856 0.34402 0 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 1339 on 4962 degrees of freedom
##
## AIC: 1355 BIC: 1407 (Smaller is better.)
#let's try a transitivity term...
# mod.tran<-ergm(g~edges+nodematch("sex")+nodematch("grade")
# +nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
# +mutual+triangle)
# summary(mod.gendmut)
#turns out that won't run...lets cheat!
mod.tran<-ergm(g~edges+nodematch("sex")+nodematch("grade")
+nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
+mutual+triangle, control=control.ergm(MCMLE.maxit=2))
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 2:
## Optimizing with step length 0.312988763891207.
## The log-likelihood improved by 5.029.
## Iteration 2 of at most 2:
## Optimizing with step length 0.000265821143928091.
## The log-likelihood improved by 0.2442.
## MCMLE estimation did not converge after 2 iterations. The estimated coefficients may not be accurate. Estimation may be resumed by passing the coefficients as initial values; see 'init' under ?control.ergm for details.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
summary(mod.tran)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") +
## nodeicov("indegree") + nodeocov("outdegree") + mutual + triangle
##
## Iterations: 2 out of 2
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -8.025142 0.405900 14 < 1e-04 ***
## nodematch.sex 0.475142 0.487653 2 0.329933
## nodematch.grade 2.195880 0.634579 2 0.000544 ***
## nodematch.smokes 0.852256 0.583772 3 0.144378
## nodeicov.indegree 0.211418 0.182808 1 0.247531
## nodeocov.outdegree 0.248641 0.162040 1 0.124985
## mutual 1.557091 0.344459 15 < 1e-04 ***
## triangle 0.106772 0.001898 7 < 1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 1468 on 4962 degrees of freedom
##
## AIC: 1484 BIC: 1536 (Smaller is better.)
mcmc.diagnostics(mod.tran)
## Sample statistics summary:
##
## Iterations = 16384:1063936
## Thinning interval = 1024
## Number of chains = 1
## Sample size per chain = 1024
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## edges 4219.3 58.02 1.8130 56.161
## nodematch.sex 2065.4 31.12 0.9726 29.214
## nodematch.grade 613.5 10.06 0.3143 8.183
## nodematch.smokes 2870.7 46.91 1.4659 44.743
## nodeicov.indegree 18093.8 128.82 4.0257 112.978
## nodeocov.outdegree 18390.6 130.86 4.0894 112.706
## mutual 2175.4 29.10 0.9094 30.775
## triangle 395589.9 7709.43 240.9196 7509.066
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## edges 4118 4132.5 4253 4255.0 4258
## nodematch.sex 2010 2023.2 2083 2084.0 2088
## nodematch.grade 595 606.5 619 619.2 621
## nodematch.smokes 2789 2801.2 2898 2899.0 2903
## nodeicov.indegree 17861 17919.2 18165 18173.0 18190
## nodeocov.outdegree 18155 18194.0 18463 18471.0 18490
## mutual 2126 2129.0 2193 2193.0 2194
## triangle 382591 382775.5 400279 400281.0 400299
##
##
## Sample statistics cross-correlations:
## edges nodematch.sex nodematch.grade
## edges 1.0000000 0.9992569 0.9932355
## nodematch.sex 0.9992569 1.0000000 0.9945745
## nodematch.grade 0.9932355 0.9945745 1.0000000
## nodematch.smokes 0.9998830 0.9992805 0.9935824
## nodeicov.indegree 0.9988885 0.9985409 0.9934550
## nodeocov.outdegree 0.9989202 0.9982089 0.9909539
## mutual 0.9978185 0.9959606 0.9883397
## triangle 0.9974954 0.9953521 0.9866779
## nodematch.smokes nodeicov.indegree nodeocov.outdegree
## edges 0.9998830 0.9988885 0.9989202
## nodematch.sex 0.9992805 0.9985409 0.9982089
## nodematch.grade 0.9935824 0.9934550 0.9909539
## nodematch.smokes 1.0000000 0.9987943 0.9987190
## nodeicov.indegree 0.9987943 1.0000000 0.9975310
## nodeocov.outdegree 0.9987190 0.9975310 1.0000000
## mutual 0.9971929 0.9950552 0.9966797
## triangle 0.9968425 0.9944679 0.9963066
## mutual triangle
## edges 0.9978185 0.9974954
## nodematch.sex 0.9959606 0.9953521
## nodematch.grade 0.9883397 0.9866779
## nodematch.smokes 0.9971929 0.9968425
## nodeicov.indegree 0.9950552 0.9944679
## nodeocov.outdegree 0.9966797 0.9963066
## mutual 1.0000000 0.9998440
## triangle 0.9998440 1.0000000
##
## Sample statistics auto-correlation:
## Chain 1
## edges nodematch.sex nodematch.grade nodematch.smokes
## Lag 0 1.0000000 1.0000000 1.0000000 1.0000000
## Lag 1024 0.9981101 0.9977836 0.9970511 0.9980660
## Lag 2048 0.9960372 0.9954138 0.9941653 0.9959323
## Lag 3072 0.9938211 0.9929642 0.9914632 0.9936546
## Lag 4096 0.9914587 0.9904576 0.9886689 0.9912237
## Lag 5120 0.9889850 0.9879250 0.9857973 0.9886773
## nodeicov.indegree nodeocov.outdegree mutual triangle
## Lag 0 1.0000000 1.0000000 1.0000000 1.0000000
## Lag 1024 0.9974614 0.9973679 0.9982533 0.9983115
## Lag 2048 0.9948052 0.9948389 0.9963786 0.9964605
## Lag 3072 0.9921616 0.9923010 0.9943711 0.9944328
## Lag 4096 0.9894583 0.9896352 0.9922170 0.9922461
## Lag 5120 0.9866739 0.9869925 0.9899313 0.9899126
##
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## edges nodematch.sex nodematch.grade
## -185.35 -112.50 -33.57
## nodematch.smokes nodeicov.indegree nodeocov.outdegree
## -245.53 -117.93 -99.36
## mutual triangle
## -279.57 -7890.08
##
## Individual P-values (lower = worse):
## edges nodematch.sex nodematch.grade
## 0.000000e+00 0.000000e+00 4.930723e-247
## nodematch.smokes nodeicov.indegree nodeocov.outdegree
## 0.000000e+00 0.000000e+00 0.000000e+00
## mutual triangle
## 0.000000e+00 0.000000e+00
## Joint P-value (lower = worse): 1.491166e-87 .
## Package latticeExtra is not installed. Falling back on coda's default MCMC diagnostic plots.
##
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).
test=simulate.ergm(mod.tran,1)
plot(test)
mod.gwesp<-ergm(g~edges+nodematch("sex")+nodematch("grade")
+nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
+mutual+gwesp)
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20:
## Optimizing with step length 0.564907830439016.
## The log-likelihood improved by 3.669.
## Iteration 2 of at most 20:
## Optimizing with step length 0.952429948637443.
## The log-likelihood improved by 2.673.
## Iteration 3 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.1513.
## Step length converged once. Increasing MCMC sample size.
## Iteration 4 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.05256.
## Step length converged twice. Stopping.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
summary(mod.gwesp)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") +
## nodeicov("indegree") + nodeocov("outdegree") + mutual + gwesp
##
## Iterations: 4 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -8.02670 0.39054 0 <1e-04 ***
## nodematch.sex 0.50762 0.11943 0 <1e-04 ***
## nodematch.grade 2.20383 0.17854 0 <1e-04 ***
## nodematch.smokes 0.78078 0.14779 0 <1e-04 ***
## nodeicov.indegree 0.23866 0.02982 0 <1e-04 ***
## nodeocov.outdegree 0.27343 0.03474 0 <1e-04 ***
## mutual 2.05411 0.24294 0 <1e-04 ***
## gwesp 0.64075 0.14652 0 <1e-04 ***
## gwesp.decay 0.06569 0.15439 0 0.671
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 6890 on 4970 degrees of freedom
## Residual Deviance: 1297 on 4961 degrees of freedom
##
## AIC: 1315 BIC: 1373 (Smaller is better.)
mcmc.diagnostics(mod.gwesp)
## Sample statistics summary:
##
## Iterations = 16384:4209664
## Thinning interval = 1024
## Number of chains = 1
## Sample size per chain = 4096
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## edges -2.1826 18.993 0.2968 1.0084
## nodematch.sex -1.2117 14.086 0.2201 0.7904
## nodematch.grade -2.1377 14.659 0.2290 0.8521
## nodematch.smokes -0.9524 15.984 0.2497 0.8639
## nodeicov.indegree -9.9270 110.148 1.7211 5.7165
## nodeocov.outdegree -19.4607 108.374 1.6933 5.6268
## mutual -1.3481 7.733 0.1208 0.4565
## gwesp -3.3380 22.683 0.3544 1.2833
## gwesp.decay -1.6761 13.610 0.2127 0.7583
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## edges -39.00 -15.00 -2.000 11.000 35.00
## nodematch.sex -29.00 -11.00 -1.000 9.000 26.00
## nodematch.grade -30.00 -12.00 -2.000 8.000 27.00
## nodematch.smokes -31.00 -12.00 -1.000 10.000 31.00
## nodeicov.indegree -225.00 -87.00 -9.000 65.000 210.00
## nodeocov.outdegree -230.62 -91.00 -20.000 54.000 197.00
## mutual -17.00 -6.00 -1.000 4.000 14.00
## gwesp -47.02 -18.36 -3.720 11.907 43.07
## gwesp.decay -27.20 -10.88 -2.058 7.213 26.51
##
##
## Sample statistics cross-correlations:
## edges nodematch.sex nodematch.grade
## edges 1.0000000 0.7920421 0.8285412
## nodematch.sex 0.7920421 1.0000000 0.6338130
## nodematch.grade 0.8285412 0.6338130 1.0000000
## nodematch.smokes 0.8821243 0.6993325 0.7473641
## nodeicov.indegree 0.9243292 0.7224975 0.6880821
## nodeocov.outdegree 0.9485801 0.7264189 0.7353055
## mutual 0.7850074 0.6248356 0.7928949
## gwesp 0.9225984 0.7214260 0.8427424
## gwesp.decay 0.7798153 0.5812022 0.7877265
## nodematch.smokes nodeicov.indegree nodeocov.outdegree
## edges 0.8821243 0.9243292 0.9485801
## nodematch.sex 0.6993325 0.7224975 0.7264189
## nodematch.grade 0.7473641 0.6880821 0.7353055
## nodematch.smokes 1.0000000 0.7572631 0.8029944
## nodeicov.indegree 0.7572631 1.0000000 0.8845197
## nodeocov.outdegree 0.8029944 0.8845197 1.0000000
## mutual 0.7077234 0.7101292 0.7345627
## gwesp 0.8103180 0.8669983 0.8927071
## gwesp.decay 0.6894998 0.7446771 0.7637741
## mutual gwesp gwesp.decay
## edges 0.7850074 0.9225984 0.7798153
## nodematch.sex 0.6248356 0.7214260 0.5812022
## nodematch.grade 0.7928949 0.8427424 0.7877265
## nodematch.smokes 0.7077234 0.8103180 0.6894998
## nodeicov.indegree 0.7101292 0.8669983 0.7446771
## nodeocov.outdegree 0.7345627 0.8927071 0.7637741
## mutual 1.0000000 0.8189673 0.7645249
## gwesp 0.8189673 1.0000000 0.8506887
## gwesp.decay 0.7645249 0.8506887 1.0000000
##
## Sample statistics auto-correlation:
## Chain 1
## edges nodematch.sex nodematch.grade nodematch.smokes
## Lag 0 1.0000000 1.0000000 1.0000000 1.0000000
## Lag 1024 0.7356949 0.7243649 0.8194960 0.7344987
## Lag 2048 0.5964351 0.5895078 0.7048158 0.6006431
## Lag 3072 0.5078308 0.4968075 0.6115425 0.5079153
## Lag 4096 0.4374859 0.4333647 0.5375935 0.4377303
## Lag 5120 0.3812865 0.3798584 0.4708243 0.3895060
## nodeicov.indegree nodeocov.outdegree mutual gwesp
## Lag 0 1.0000000 1.0000000 1.0000000 1.0000000
## Lag 1024 0.7273865 0.7310531 0.8378584 0.7870712
## Lag 2048 0.5831148 0.5891169 0.7203466 0.6655837
## Lag 3072 0.4913930 0.4989458 0.6233013 0.5797273
## Lag 4096 0.4242818 0.4363518 0.5490370 0.5117562
## Lag 5120 0.3672405 0.3829011 0.4868593 0.4542104
## gwesp.decay
## Lag 0 1.0000000
## Lag 1024 0.7811626
## Lag 2048 0.6533085
## Lag 3072 0.5649845
## Lag 4096 0.4998976
## Lag 5120 0.4413901
##
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## edges nodematch.sex nodematch.grade
## 0.5028 1.3126 1.1333
## nodematch.smokes nodeicov.indegree nodeocov.outdegree
## 0.2698 0.2840 0.3324
## mutual gwesp gwesp.decay
## 1.5353 0.8789 0.8636
##
## Individual P-values (lower = worse):
## edges nodematch.sex nodematch.grade
## 0.6151170 0.1893017 0.2570862
## nodematch.smokes nodeicov.indegree nodeocov.outdegree
## 0.7872953 0.7764182 0.7396138
## mutual gwesp gwesp.decay
## 0.1246995 0.3794302 0.3877901
## Joint P-value (lower = worse): 0.09919254 .
## Package latticeExtra is not installed. Falling back on coda's default MCMC diagnostic plots.
##
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).
test=simulate.ergm(mod.gwesp,1)
plot(test)
try(plot(gof(mod.gwesp), main=NULL))
################################################################
#latent space models. These models allow us to fit a network without
#having to specify the actual tie formation patterns...
#code cribbed from Jake Fisher
#################################################################
# The most user-friendly latent space model software is in the latentnet package
# (more recent models are provided by the amen package).
#install.packages("latentnet")
library(latentnet)
##
## latentnet: version 2.8.0, created on 2017-08-20
## Copyright (c) 2017, Pavel N. Krivitsky, University of Wollongong
## Mark S. Handcock, University of California -- Los Angeles
## with contributions from
## Susan M. Shortreed
## Jeremy Tantrum
## Peter D. Hoff
## Li Wang
## Kirk Li, University of Washington
## Jake Fisher
## Jordan T. Bates
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("latentnet").
## NOTE: BIC calculation prior to latentnet 2.7.0 had a bug in the calculation of the effective number of parameters. See help(summary.ergmm) for details.
## NOTE: Prior to version 2.8.0, handling of fixed effects for directed networks had a bug. See help("ergmm-terms") for details.
start.time <- Sys.time()
latent.fit <- ergmm(g ~ euclidean(d = 2))
runtime=Sys.time()-start.time;
runtime;
## Time difference of 1.408969 mins
summary(latent.fit)
## NOTE: It is not certain whether it is appropriate to use latentnet's BIC to select latent space dimension, whether or not to include actor-specific random effects, and to compare clustered models with the unclustered model.
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ euclidean(d = 2)
## Attribute: edges
## Model: Bernoulli
## MCMC sample of size 4000, draws are 10 iterations apart, after burnin of 10000 iterations.
## Covariate coefficients posterior means:
## Estimate 2.5% 97.5% 2*min(Pr(>0),Pr(<0))
## (Intercept) 1.7690 1.3795 2.1624 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Overall BIC: 1910.635
## Likelihood BIC: 1127.11
## Latent space/clustering BIC: 783.5258
##
## Covariate coefficients MKL:
## Estimate
## (Intercept) 0.9692998
plot(latent.fit)
mcmc.diagnostics(latent.fit)
## Chain 1
## Lag 0
## lpY beta.1 Z.1.1 Z.1.2
## lpY 1.0000000 0.4594187 0.1319291 0.1317427
## beta.1 0.4594187 1.0000000 0.1526472 0.2620373
## Z.1.1 0.1319291 0.1526472 1.0000000 0.2264011
## Z.1.2 0.1317427 0.2620373 0.2264011 1.0000000
##
## Lag 10
## lpY beta.1 Z.1.1 Z.1.2
## lpY 0.45622566 0.29748463 0.10352657 0.07067086
## beta.1 0.30626407 0.31471132 0.06964949 0.06954152
## Z.1.1 0.09422899 0.05918816 0.51834688 0.09056140
## Z.1.2 0.07853141 0.06887424 0.10328924 0.45935555
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.0125
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## lpY 40 7800 600 13.0
## beta.1 30 7170 600 12.0
## Z.1.1 40 7970 600 13.3
## Z.1.2 100 17340 600 28.9
plot(gof(latent.fit), main=NULL)
# Can add additional dimensions...
start.time <- Sys.time()
latent.fit <- ergmm(g ~ euclidean(d = 3))
runtime=Sys.time()-start.time;
runtime;
## Time difference of 1.981795 mins
plot(gof(latent.fit), main=NULL)
# ... latent groups ...
start.time <- Sys.time()
latent.fit <- ergmm(g ~ euclidean(d = 2, G = 5))
runtime=Sys.time()-start.time;
runtime;
## Time difference of 2.776853 mins
plot(latent.fit)
plot(gof(latent.fit), main=NULL)
# ... or homophily effects
start.time <- Sys.time()
latent.fit <- ergmm(g ~ nodematch("grade") + nodematch("sex")+euclidean(d = 2))
runtime=Sys.time()-start.time;
runtime
## Time difference of 1.369521 mins
summary(latent.fit)
##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: g ~ nodematch("grade") + nodematch("sex") + euclidean(d = 2)
## Attribute: edges
## Model: Bernoulli
## MCMC sample of size 4000, draws are 10 iterations apart, after burnin of 10000 iterations.
## Covariate coefficients posterior means:
## Estimate 2.5% 97.5% 2*min(Pr(>0),Pr(<0))
## (Intercept) -1.11165 -1.58619 -0.6259 <2e-16 ***
## nodematch.grade 3.67247 3.23129 4.0931 <2e-16 ***
## nodematch.sex 0.27219 -0.20578 0.6836 0.223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Overall BIC: 1741.307
## Likelihood BIC: 1161.369
## Latent space/clustering BIC: 579.9382
##
## Covariate coefficients MKL:
## Estimate
## (Intercept) -1.9357062
## nodematch.grade 3.2004359
## nodematch.sex 0.2859632